/*
    Copyright 2006-2012 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains classes and functions related to blackbody radiation.

#ifndef __blackbody__
#define __blackbody__

#include "mcrx-types.h"
#include "constants.h"
#include "blitz/array.h"

namespace mcrx {
  class blackbody;

  T_float B_lambda (T_float, T_float);
  T_float dB_lambda_dT (T_float, T_float);
  T_float dB_lambda_dinvT (T_float, T_float);

   BZ_DEFINE_BINARY_FUNC(Fn_blambda, B_lambda);
   BZ_DECLARE_ARRAY_ET_BINARY(B_lambda, Fn_blambda);
   BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(B_lambda, Fn_blambda, T_float);

   BZ_DEFINE_BINARY_FUNC(Fn_dblambda, dB_lambda_dT);
   BZ_DECLARE_ARRAY_ET_BINARY(dB_lambda_dT, Fn_dblambda);
   BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(dB_lambda_dT, Fn_dblambda, T_float);

   BZ_DEFINE_BINARY_FUNC(Fn_dbilambda, dB_lambda_dinvT);
   BZ_DECLARE_ARRAY_ET_BINARY(dB_lambda_dinvT, Fn_dbilambda);
   BZ_DECLARE_ARRAY_ET_BINARY_SCALAR(dB_lambda_dinvT, Fn_dbilambda, T_float);
}

/** Calculates the Planck function for a blackbody of temperature T at
    wavelength 1/invlambda. The inverse wavelength and temperature is
    passed for efficiency. */
inline mcrx::T_float mcrx::B_lambda (T_float invlambda, T_float invT) 
{
  assert(invlambda>0);
  assert(invT>0);

  // 2hc^2 in SI units [kg m^4 s^-3].
  const T_float a = 1.19104393407e-16;
  // hc/k for the exp function, in SI units [m K].
  const T_float d = 1.43876866033e-2;
  const T_float il2 = invlambda*invlambda;
  const T_float e = expm1 (d*invlambda*invT);

  const T_float ret=a*il2*il2*invlambda/e;
  assert(ret>=0);
  return ret;
}

/** Calculates the temperature derivative of the Planck
    function for a blackbody of temperature 1/invT at wavelength
    1/invlambda. The inverse wavelength and temperature is passed for
    efficiency. */
inline mcrx::T_float mcrx::dB_lambda_dT (T_float invlambda, T_float invT) 
{
  // This is almost identical do the inverse-temp derivative, so we
  // just use that function.
  return -dB_lambda_dinvT(invlambda, invT)*invT*invT;
}


/** Calculates the inverse-temperature derivative of the Planck
    function for a blackbody of temperature 1/invT at wavelength
    1/invlambda. The inverse wavelength and temperature is passed for
    efficiency. */
inline mcrx::T_float mcrx::dB_lambda_dinvT (T_float invlambda, T_float invT) 
{
  assert(invlambda>0);
  assert(invT>0);

  // 2hc^2 in SI units [kg m^4 s^-3].
  const T_float a = 1.19104393407e-16;
  // hc/k for the exp function, in SI units [m K].
  const T_float d = 1.43876866033e-2;
  // save the exponent expression
  //const T_float e = exp (d*invlambda*invT);
  const T_float iem1 = 1./expm1 (d*invlambda*invT);
  const T_float il2 = invlambda*invlambda;

  // we should have e/em1^2, but it is ill behaved when e is
  // inf. Instead we use (1/em1 + 1/em1^2) which is the same thing.
  const T_float ret = -a*d*il2*il2*il2*(iem1+iem1*iem1);
  assert(ret==ret);
  return ret;
}


/** Calculates the emission from a blackbody with a specified
    temperature and luminosity. Can be used to set luminosity of the
    emitters to a BB, but this by itself is just a utility class. */
class mcrx::blackbody {
private:
  T_float temperature_; ///< The temperature of the black body.
  T_float luminosity_; ///< The luminosity of the BB.

public:
  blackbody(T_float t, T_float l): temperature_ (t), luminosity_ (l) {};

  T_float emission(T_float lambda) const;
  array_1 emission(const array_1& lambda) const;

  T_float temperature () const {return temperature_;};
  T_float luminosity () const {return luminosity_;};
};

#endif
